Thermomechanical properties of graphene: valence force field model approach 
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Using the valence force field model of Perebeinos and Tersoff [Phys. Rev. B 79, 241409(R) (2009)] , 
different energy modes of suspended graphene subjected to tensile or compressive strain are studied. 
By carrying out Monte Carlo simulations it is found that: i) only for small strains (|e| ^ 0.02) 
the total energy is symmetrical in the strain, while it behaves completely different beyond this 
threshold; ii) the important energy contributions in stretching experiments are stretching, angle 
bending, out-of-plane term and a term that provides repulsion against tt — tt misalignment; 
iii) in compressing experiments the two latter terms increase rapidly and beyond the buckling 
transition stretching and bending energies are found to be constant; iv) from stretching-compressing 
simulations we calculated the Young modulus at room temperature 350±3.15 N/m, which is in good 
agreement with experimental results (340±50 N/m) and with ab-initio results [322-353] N/m; v) 
molar heat capacity is estimated to be 24.64 J/mol~^K~^ which is comparable with the Dulong- 
Petit value, i.e. 24.94 J/mol~^K~^ and is almost independent of the strain; vi) non-linear scaling 
properties are obtained from height-height correlations at finite temperature; vii) the used valence 
force field model results in a temperature independent bending modulus for graphene, and viii) the 
Gruneisen parameter is estimated to be 0.64. 

Keywords : Thermomecahnical properties. Strain graphene sheet, Monte Carlo simulation. Molar 
heat capacity 
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I. INTRODUCTION 

Since the discovery of graphene in 2004, which is an al- 
most two dimensional crystalhne material, its exceptional 
mechanical properties have been studied [FJiJ. Ten- 
sional strain in monolayer graphene affects its electronic 
structure. For example strains larger than 15% changes 
graphene's band structure and leads to the opening of 
an electronic gap [7[. In recent experiments the buck- 
ling strain of a graphene sheet that was positioned on 
top of a substrate was found to be six orders of magni- 
tude larger (i.e. 0.5 — 0.6%) than for graphene suspended 
in air j5i] . Furthermore, some experiments showed that a 
compressed rectangular monolayer of graphene on a plas- 
tic beam with size 30x100 ^m^ is buckled at about 0.7% 
strain [81]. 

Elasticity theory for a thin continuum plate and the 
empirical interatomic potentials (EP) are two main the- 
oretical approaches that have been used to study vari- 
ous mechanical properties measured in compressing and 
stretching experiments [J, M, ll3- Continuum elas- 
ticity theory does not give the atomistic features of 
graphene while the EPs, such as the Brenner poten- 
tial (REBO) [nl, [13 and the LCBOPII potential [H, 
can properly account for the mechanical properties of 
graphene. Despite the several benefits of these EPs, some 
special atomistic features of graphene subjected to com- 
pressive or tensile strains could not be explained. The dif- 
ferent energy contributions in these potentials are mixed. 
For example, in REBO all the many body effects are put 
in the bond order term and the different important en- 
ergy contributions are not separable. 

It still remains unclear how large are the contributions 



of the different energy terms in strained graphene. Using 
the recently introduced valence force field (VFF) model 
by Perebeinos and Tersoff jlj] we show how the contribu- 
tion of the different energy modes in strained graphene 
can be separated and we calculate their dependence on 
the value of strain. 

The bending modulus of graphene at zero tempera- 
ture was estimated using several interatomic potentials, 
e.g. the first version of the Brenner potential [II5 yields 
0.83 eV, the second generation of the Brenner poten- 
tial [12] estimated it to be 0.69 eV, adding the third near- 
est neighbors (the dihedral angle effect) in the Brenner 
potential enhances it to 1.4 eV [15|, using the LCBOPII 
potential and continuum membrane theory the bending 
rigidity was found to be 0.82 eV d [H, Tersoff 's VFF 
model estimated it to be 2.1 eV, and from ab-initio en- 
ergy calculations it was found to be 1.5 eV [16] (note 
that 'bending rigidity' ('bending modulus') is used for a 
membrane stiffness (an atomistic sheet)). Despite these 
studies the temperature dependence of the bending mod- 
ulus is poorly known. An increasing behavior versus 
temperature for the bending rigidity was found [4] by 
using Monte Carlo simulations with the LCBOPII poten- 
tial and membrane theory concepts. In contrast, Liu et 
al [l7[ found a decreasing bending rigidity with tempera- 
ture using the REBO. Here we show that the VFF model 
predicts a temperature independent bending modulus. 

In this study we employ VFF and carry out standard 
Monte Carlo simulations in order to calculate and com- 
pare the different energy modes of a graphene sheet that 
is subject to axial strains. The total energy is found to 
be different for compressing and stretching when strains 
are applied larger than |2|%. Two important terms, i.e. 
stretching and bending, vary differently depending on the 



way that one stretches or compresses the system. We find 
that out-of-plane and tt — tt terms have much larger con- 
tributions in compression experiments when compared 
to stretching. Furthermore, we used this potential to 
calculate Young's modulus at room temperature from 
stretching-compressing simulations. We also calculate 
the molar heat capacity. Our Monte Carlo simulations 
show that the VFF potential yields a temperature inde- 
pendent bending modulus. 

This paper is organized as follows. Section II contains 
the essentials of the VFF model for graphene. The sim- 
ulation method for strained graphene will be presented 
in Sec. III. Different energy modes of strained graphene 
are studied in Sec. IV. In Sec. V the molar heat capacity 
for non-strained suspended graphene is calculated. Tem- 
perature effects of the bending modulus of graphene with 
periodic boundary condition are presented in Sec. VI and 
the scaling properties of graphene at finite temperature 
are investigated in Sec VII. We will conclude the paper 
in Sec. VIII. 



II. ELASTIC ENERGY OF GRAPHENE 

There are two main classical approaches for the investi- 
gation of the elastic energy of graphene: 1) the continuum 
approach based on elasticity theory, and 2) the atomistic 
description using accurate interatomic potentials. 

The total energy of a deformed membrane con- 
sists of two important terms: stretching and 
bending. For almost flat and continuum mem- 
brane using Monge representation the surface 
area element dA can be approximated by a 
flat sheet area element in the x-y-plane, i.e. 
dA « dxdy and the bending energy is written 
as ^ J dxdyK{\/^h)^ where k is the bending rigid- 
ity and h is the out-of-plane deformation of the 
membrane at point {x,y). The stretching term 
for an isotropic continuum material in the lin- 
ear regime includes two independent parameters: 
the shear modulus (/i) and the Lame coeflScient 
(A) and is written as ^ J dxdy[2fj,u'^p -I- Au^^]. Here 

Uai3 = ^[daUjj + dpUa + dahdjjh] is the second rank 
symmetric tensor with a,/? = 1,2 and Ua{x,y) 
is the a*'' component of the displacement vec- 
tor. Neglecting the last term in the strain tensor 
makes the stretching term linear and decouples 
the bending and stretching energy. Therefore for 
an isotropic and continuum material for small de- 
formations and with the assumption of a nearly 
flat membrane (|V/ip ^ 1) the strain energy (Ut) 
can be written as [18] 



1 



Ut ^ - dxdy[K{\I^hY + 2^w„^ + At 



(1) 



The integral is taken over the projected area of the mem- 
brane into the x-y-plane. For isotropic materials and in 
the linear approximation the mentioned parameters are 



related to the Young modulus (Y) and Poisson's ratio {v) 
as /i = 1"/(2(1 -f v)) and A = 2/ij//(l - 2v). Equation (P) 
can be rewritten in terms of the Fourier components of 
h and yields the scaling properties of the sheet. Despite 
these benefits, this continuum model does not include 
self-avoidance, the natural condition of true physical sys- 
tems and does not show atomistic details of the mem- 
brane under different boundary conditions. All these de- 
ficiencies originate from the continuity assumption. As- 
suming graphene as a continuum plate limits the study 
to only bending and stretching modes. 

Due to the hexagonal symmetry of the flat 
monolayer graphene lattice, it is elastically 
isotropic which implies that the the bending mod- 
ulus is independent of the direction at least within 
the linear elastic regime [15]. However, the 
graphene monolayer can exhibits anisotropic be- 
havior in the nonlinear regime where distortions 
are no longer infinitesimal. The larger stretches, 
the stronger anisotropy and non-linearity eff'ects. 
Cadelano et al found that monolayer graphene 
is isotropic in the linear regime, while it is 
anisotropic when nonlinear features are taken 
into account [19.] . 

The recently introduced VFF model in Ref. [Til] is 
expected to be able to describe both compression and 
stretching experiments by separating the contribution of 
the various energy modes. This model includes explic- 
itly the various relevant energy terms which describe the 
change in the bond lengths, bond angles and torsional 
effects. The total energy density is written as 
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NSi 



■{Est + Ebe + Eout + Ebo + Ep + Eco), (2) 



where N is the number of atoms and Sq — ^^a§ is the 
surface area of the unit cell of the honeycomb lattice. In 
the following we will discuss the different terms in Eq. 
(2). Note that the energy reference is set to zero. As- 
suming flo — 1.42 A as the unit of length, the ^stretchings 
and ^bending' (bending of the bond angle) terms are 



Est = \KsY,{5n,f 



(3) 



Ebe = Kie J2 (C0s(%fc) - COS(0O))' (4) 

where Svij — Vij — 1 and ^o — 27r/3. In Eqs. (3) and 
(4) rij is the bond length between atom 'i' and 'j', 9ijk 
is the angle between the nearest neighbor atoms 'i', 'j' 
and 'k' and 9o is the equilibrium angle between three 
nearest neighbor atoms. Est is the two-body stretching 
term responsible for bond stretching. Ei,e is the bending 
energy due to the bond angles. Here, all bond angles will 
be considered. The above two terms results in a quasi- 
harmonic model |l3i. Later, we will find that these two 



TABLE I: Parameters of the energy model (Eq. ([2])) are taken 
from Ref. 



M- 


Jnits 


ire in 


eV. 






Ks 


Kbe 


Kout 


Kbo 


Kr, 


Kcor 


37.04 


4.087 


1.313 


4.004 


0.016102 


4.581 



terms become constant as function of strain when beyond 
the buckhng point in a compression experiment. 

The stiffness against ^ out-of-plane^ vibration is pro- 
vided by 



Env.t — K„ 



E<; 



'^'^ij ' '^ik ^ '^il 






)^ (5) 



where the summation is taken over the first neighbors 
of atoms 'z' and taking care of not double counting. In 
Eq. (5) rij is the distance vector between atom 'i' and 
'j'. Hence there are three different terms for each atom. 
Correlations between bond lengths are provided by the 
^bond order' term 



Ebo = Kbo ^ Sr^jSr^k, 

i^j<k 



(6) 



where for each bond length with central atom 'i' three 
different terms are considered. 

The misalignment of the neighboring n orbital is given 
by the 'tt — tt' term 



Ep ^ l^Kp^l'^^ X TT 



^A\ 



i,3 



where 



^ = 3- 






(7) 



(8) 



nijk — Tij X rik is a vector normal to the plane passing 
through the vectors rij and rik- This kind of interac- 
tion plays an important role in the interlayer interaction 
in graphitic structures. Note that the simple two body 
interaction gives only 2% of the local density of state 
(LDA) result for the energy difference between AA and 
AB stacked graphite [20| . 

The last term takes into account the coupling between 
bond stretching and bond angle bending (bond length- 
bond angle cross coupling), i.e. the 'correlation' term 

Ecor = Kcor ^ (5r.y (cos(%fe) - cos(6lo)). (9) 

The coefficients in the above equations {Ks,Kbe and 
so on) were recently parameterized by Perebeinos and 
Tersoff [ij] such that the phonon dispersion of graphenc 
was accurately described. These parameters are listed in 
Table I. 



III. SIMULATION METHOD: 
GRAPHENE 



STRAINED 



In order to compress (stretch) graphene nanoribbons, 
we have carried out several standard Monte Carlo simula- 
tions [J| of a suspended graphene sheet at finite tempera- 
ture. Equation ((2]) is used to calculate the total energy of 
the system. Our sample is a rectangular graphene sheet 
with Ix X ly dimensions in x- and y-directions, respec- 
tively, containing 7Vo=1600 atoms. The sheet is strained 
along the armchair or the zig-zag direction. Strain is al- 
ways applied along x. When graphene is strained along 
the armchair direction we named it armchair graphene 
-AC- (^a;=85.2A, ly=4Q.l%k) and strained graphene 
along zig-zag direction is named zig-zag graphene-ZZ- 
(Zj;=98.38A, Zy=42.6A). Periodic boundary conditions 
are applied along the lateral direction i.e. zigzag direc- 
tion in AC graphene and along the armchair direction 
for ZZ graphene. Our simulation starts with a flat sheet, 
and we allow then the system to thermally equilibrate 
such that the total energy no longer changes. Tempera- 
ture is typically taken T=300K, except when otherwise 
indicated. Figure [ija) shows a snap shot of the relaxed 
unstrained ZZ graphene at T=300K (note that the sup- 
ported ends are fixed) . We found that the graphene sheet 
is corrugated after relaxation which are the intrinsic ther- 
mal ripples in graphene. Thus the used VFF is able to 
display true structural properties. These ripples are vi- 
tal in order to make suspended graphene stable and are 
therefore crucial for the stability of a fiat 2D crystal at 
finite temperature [J]. 

To simulate a suspended sheet we fixed two atomic 
rows at both longitudinal ends. These boundary atoms 
are not included in the summations when calculating 
the different energy contributions (in Eqs. (2-9)), i.e. 
N — Nq — 40. We compress/stretch the system with 
a slow rate, i.e. in every million Monte Carlo steps 
the longitudinal ends are reduced/elongated with about 
5 =0.02 A such that the system stays in thermal equilib- 
rium. After obtaining the total desired strain e, we wait 
for an extra 4 million steps during which the system can 
relax. For example, a strain (applied in x-direction) of 
e = 1.2% is achieved after 29x10^ Monte Carlo steps. 



IV. DIFFERENT ENERGY MODES FOR 
STRAINED GRAPHENE 

Figures [ljb,c) show two snap shots of compressed ZZ 
and AC graphene, respectively when e=-2%. It is inter- 
esting that the rippled structure is different in the two 
cases. This is due to the different out-of-plane and tt — tt 
interaction terms around and beyond the buckling tran- 
sition points, i.e. e < —2.5%. 

The variation of height, Ah = a/< h'^ > - < h >2, 
in Fig. 1(a) after 10 million MC steps fluctu- 
ates around 0.2 A which is comparable with those 
found when using REBO [21.] . In Figs. l(b,c) for 





TABLE II: Young's modulus of graphene in units of N/m. 




Experimental 


Classical (T=300 K) 


Ab-initio (T=OK) 


Tight-Binding 


340±50 


350±3.15 355±21 384 


345±6.9 336 352.54 351.75 322 


312 


Ref. yj 


Present* Ref. flO]" Ref. [22]"^ 


Ref. fl6] Ref. [23] Ref. [24] Ref. [25] Ref. [26] 


Ref. [19] 



* VFF model [ij, '' LCBOPII potential, " Tersoff-Brenner potential 




FIG. 1: (Color online) Snap shot of a suspended graphene 
sheet at T=300K using the valence force model (Eq. (djl). 
Blue lines indicate the position of fixed atoms in a; — y plane, 
(a) Unstrained, (b) compressed ZZ graphene, and (c) com- 
pressed AC graphene. 



compressed nanoribbons of about e = -2.0 % A/i 
is 0.5 A after 54 million MC steps. The larger 
compressive strain yields a larger height variance. 

Figure [2Ja) shows the variation of the total energy 
(Eq. ^) with applied strain at T=300K. The verti- 
cal dashed line separates compressive (left) and tensile 
strain (right). Square (circular) symbols refer to AC (ZZ) 
graphene. Notice that AC and ZZ strained graphene re- 
sult in the same energy, although their ripples structure 
(see Figs. [llb,c)) can be rather different. Note that the 
energy curve is no longer symmetric around e = beyond 
the colored rectangle where |e| > 0.02. Inside this region 
the deformation is symmetric and the harmonic approx- 
imation to the total energy works well as shown by the 
full black (parabolic) curve in Fig. [2ja). The solid curve 
is a quadratic fit according to Et = Eq + \Ye'^ for only 
positive strains, where the fitting parameter Y is Young's 
modulus and Eq is the energy of the graphene sheet in 



the absence of strain. We found £;o=0.232±0.002N/m 
and Y=350.42±3.15N/m for room temperature. The 
calculated error bars are derived from the fit- 
ting procedure of our numerical data. The best 
fit yielded the smallest deviation from the har- 
monic behavior. Our result for the room tempera- 
ture Young's modulus is close to the experimental value 
(340±50N/m) and is within the ab-initio results (335- 
353 N/m) and is in agreement with those obtained 
from other classical force fields such as (LCBOPII [10l | 
and Tersoff-Brenner [24I) and Tight-Binding [13|, see 
Table. II. Note that Perebeinos and Tersoff estimated 
Y at zero temperature and found 1.024 N/m^ (343.04 
N/m) [ij]. Here we calculated Y at room temperature 
via stretching-compression simulations. Different force 
fields are parameterized such that they describe a 
set of chosen experimental data of particular ex- 
perimental effects. For example, the VFF model 
can not be used to study hydrogenation, melting 
and defect formation in either graphene or car- 
bon nanotubes sheets, while the REBO has been 
set-up such that it can be used in those cases. 
The property that the energy can be separated 
into different energy modes and the simplicity of 
coding the VFF potential are two important ad- 
vantages of this model. 

Notice that the total energy for AC (square symbols 
in Fig. [2Ia)) and ZZ (circular symbols in Fig. [UJa)) 
graphene are almost the same which is in agreement with 
the results of Ref [l^ . Graphene acts isotropically in the 
linear elastic limit. Beyond the harmonic regime there is 
a small local maximum in the energy for compression 
which is related to the buckling of graphene. Notice that 
in this regime there are small differences between ZZ and 
AC sheets. The buckling threshold is about e^ ~ —2.5%. 
The buckling strains is smaller than those found by us- 
ing REBO 19], i.e. -0.86%. Notice that both the bound- 
ary conditions and the employed interatomic potentials 
are responsible for the difference in the buckling thresh- 
olds. The main difference is due to the different potential. 
The VFF model is not a bond-order potential (REBO). 
As we mentioned in the introduction the bending mod- 
ulus predicted by REBO is about 0.69-0.83 eV which is 
smaller than the one predicted by the VFF model (2.1 
eV). Therefore, we expect a larger buckling transition us- 
ing REBO and a smaller one using VFF model (consid- 
ering the negative sign for compressive strains). 
Another important reason for the different result is the 
calculations method. Here we used Monte Carlo (time is 
meaningless) and in Ref. [9] we used Molecular dynamics 
simulations. 
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FIG. 2: (Color online) (a) Total energy of a graphene sheet subjected to stretching and compression for AC and ZZ. (b) 
Contribution of the bending (Eq. (4)) and the stretching (Eq. (3)) terms of the total energy for AC. (c) Contribution of the 
other remaining terms given by Eqs. (5)-(9) for AC. 



S(*i - 

4(f/a 

2 cm 
I cm 

cm 



-1CK- 




AC 



if =0.0352 
■ t-0.0164 



Kt-0.0 



Stretch 



Bending Bond order Oiil-ul-plaiie 





ZZ 



HE = 0.0366 
■ [■ = 0.0142 



<e: = o.o 



II j| J 

Slrclch Deiidiiig Dontiorder Out-of-plmie ^pP^ '^ 



FIG. 3: (Color online) Contribution of the different energy 
terms to the total energy for three typical values of the strain 
in AC (Top) and ZZ (Bottom) graphene. 



Figure Hlb) shows the contribution of the two impor- 
tant energy terms, i.e. stretching and bending as given 
by Eq. (3) and Eq. (4), respectively, for strained AC 
graphene. Notice that the stretching energy is larger than 
the bending and that the rate of increase for stretching is 
different. In the compression part (i.e the region to the 
left of the vertical dashed line) , after the buckling points 
these energies are almost constant. Thus increasing com- 
pression beyond the buckling point does not change the 
bending and stretching energies. Figure [3] shows the con- 
tribution of the different energy terms (scaled by Ex) for 



three values of the strain for both AC and ZZ graphene 
(e.g. the first set of bars to the left refer to 100 x Egt/ET 
for each particular strain shown in the legends). From 
Fig. [3l we conclude that the contribution of the en- 
ergy terms (Eq. (4) and Eq. (7)) which are not de- 
terminable in the continutim elasticity energy approach 
(Eq. ([Ij) are substantial and should be retained when de- 
scribing strained graphene at the atomistic scales 

Figure [SJc) shows the variation of the other terms in 
the energy, Eqs. (5-9), with strain. The energy for tt-tt 
repulsion and out-of-plane increase (decrease) with com- 
pression (stretching) and they behave opposite to the 
other terms. In compression experiments the sheets be- 
come strongly corrugated and neighbor n orbitals become 
more misaligned. In other words the normal to the ad- 
jacent surfaces, e.g. - n^fe- and -riiji- become more mis- 
aligned which results in an increase of the total energy. 
Notice also that the bond order term is smaller and neg- 
ative in the compression part with respect to the stretch- 
ing part. The correlation between the bond lengths is 
always negative for the compression part. We found that 
the relative contribution of the different energy modes 
for stretching and compression of AC and ZZ graphene 
are almost the same, compare Fig. [HJa) and Fig. ^h). 
However as we see from Figs. [Hb,c) the structure of the 
ripples depends on the direction of the applied strain. In 
the case of AC the ripples are regular (sinusoidal shape) 
while they exhibit an irregular pattern in ZZ graphene. 
Note that the buckling in plates is generally known to 
depend on the plate geometric parameters [3, |27[ . No- 
tice that, the dependence of the ripple structure of the 
graphene sheets on the sheet geometry has been demon- 
strated experimentally in Ref. Q- 



V. MOLAR HEAT CAPACITY 

Next, we simulated graphene at different tempera- 
tures. Figure Hlja) shows the temperature dependence 
of the average total energy and of the six energy terms 



for £ = 0. Notice that all energy terms vary linearly 
with T. The quantity Cy = NaSq-^^ gives the po- 
tential energy contribution to the molar heat capacity 
of the system at constant volume, where Nj^ is Avo- 
gadro's constant. Note that we first relaxed the vol- 
ume of the system by performing a constant pressure- 
temperature (NPT) Monte Carlo simulation (which re- 
moves possible boundary strains). Then we fixed the 
boundaries to the found relaxed size and allowed for 
additional thermal relaxation (i.e. constant volume- 
temperature Monte Carlo simulation or NVT). During 
this new thermal relaxation no strain is applied e = 0, 
thus, the calculated heat capacity corresponds to con- 
stant volume molar heat capacity. Surprisingly, we found 
that Cv — 12.33 Jmol~^ K~^ which is almost half of 
the Dulong-Petit value, i.e. 33fi=24.94 Jmol"^ K'^ No- 
tice that (Et) is the average of the potential energy of 
graphene which is taken over 4 million Monte Carlo steps. 
Assuming that the average of the kinetic term equals the 
average potential energy {{Et)) according to the equipar- 
tition theorem (in the harmonic regime), we can write the 
total energy {E) = {Et) + {K) = 2{Et) and then the to- 
tal heat capacity is found to be 24.66±0.10 Jmol~^ K^^. 
The obtained result is close to our previous result ob- 
tained using REBO, i.e. 24.98±0.14Jmol-iK-i Q. In 
Ref. Il0| the heat capacity at 300 K was found to be 
24.2 Jmol~^ K~^. We have performed many simulations 
at different temperatures for strained graphene and found 
always a linear {Et)—T curves. Fig.Hl^b) shows the vari- 
ation of Cy versus strain. It is interesting to note that 
Cv is slightly lower (~ 1.0%) in compressed graphene as 
compared to stretched graphene. 

Furthermore, we performed several simulations 
using the same sample employing the AIREBO 
potential [29] v^rithin LAMMPS software [30]. It 
is interesting to know that AIREBO gives (in 
the range of lOK-lOOOK) Cy=24.92 J mol ^ K ^ 
which was found to be independent of temper- 
ature. Therefore, the VFF model, REBO and 
AIREBO predicts temperature independent heat 
capacity. 

On the other hand we found that the VFF model, 
REBO and AIREBO, give a linear increase in 
the carbon-carbon bond length (a) with temper- 
ature. The resulting bond length thermal expan- 
sion coefHcients for the VFF model, REBO and 
AIREBO are a = ^^ = (5.0 ± 0.07) x IQ-^K'^ 
(5.0 ± 0.03) X 10-^K-^ and (7.0 ± 0.04) x IQ-^R-^ 
respectively. The Gruneisen parameter is de- 
fined as 7 = 7^^ where B is the two dimen- 

' Cvp 

sional bulk modulus for graphene, i.e B = 12.7 
eVA ^ [10]) s^nd p is the mass density of graphene, 
i.e. p = 12.0/5*0 = 7.6 x lO^^'gm^^. Using our re- 
sult for Cv and a gives 7 = 0.64 which is better 
estimation for the Gruneisen parameter than the 
one found in Ref. |l4l |. i.e. -0.2, and is closer to 
the experimental result, i.e. 2.0 j3lj]. 



VI. TEMPERATURE EFFECT OF THE 
BENDING MODULUS 

A common method for calculating the bending mod- 
ulus of graphene is by performing several simulations as 
function of the radius (i?) of the curved tubes and extrap- 
olates the results to i? — >■ 00 (see Fig. 5(b)). Hence, one 
can calculate the elastic energy of carbon nanotubes as a 
function of the inverse square of the radius, E — ^kR^^ . 
The coefficient k in the elastic energy gives the bend- 
ing modulus of graphene. In order to study the effect 
of temperature on the bending modulus of graphene we 
have carried out several NPT Monte Carlo simulations 
(constant pressure with periodic boundary condition) at 
different temperatures. For each particular temperature 
we have 8 different tubes. In this part of the paper, our 
systems are different armchair carbon nanotubes with ra- 
dius R — 3m^ ao/2TT and initial length 10 nm. We used 
eight armchair carbon nanotubes with index (m,m) for 
m=:5, 10, 15, 20, 25, 30, 35, 40. For each particular nanotube 
with index m, we carried out several NPT Monte Carlo 
simulations with periodic boundary condition along the 
nanotube axis and varying temperatures in the range 10 
to 1000 K. Calculating Et by using Eq. (2) for ah nan- 
otubes at temperature T, we fitted ^kR^^ to the data 
and found the bending modulus (stiffness), k, at T. From 
Fig. EKa) we notice that k is practically temperature 
independent and is about 2.02 eV. Thus the present 
VFF model results in a temperature independent bend- 
ing modulus. Using membrane theory to calculate the 
bending rigidity of graphene shows that different poten- 
tials leads to conflicting temperature dependence for the 
bending rigidity, e.g. LCBOPII [^ yields an increase 
of the bending rigidity with tem per ature while REBO 
predicts a decreasing dependence [17| . 



VII. SCALING PROPERTIES 

In the harmonic regime the power spectrum of the 
graphene solid membrane can be obtained by calculat- 
ing < |ft,qp > where hq is the Fourier transform of the 
height of the atoms {h) and q is the norm of the wave 
vector "^ (= {qx,qy) = 27r(j^, ^)) with integers Ux and 
Uy where 1^. and ly are the longitudinal and lateral size 
of the graphene sample. It is important to note that 
in this section we simulated a graphene sheet with ini- 
tial size Ix = 230.04 A and ly = 221.351 (iV=19440) 
using standard NPT Monte Carlo simulations with peri- 
odic boundary conditions in both directions (the method 
is similar to that reported in Refs. [4,10]). We esti- 
mated the spectral modes hq by fitting |ft,gp to a g" func- 
tion, from which we extract the power law a. Figure [S] 
shows the variation of \hq\'^ (averaged over 500 Monte 
Carlo realizations where 5 neighboring points were ac- 
cumulated and averaged to a single point in order to 
make the curves smoother) versus q for graphene at two 
temperatures 200 K and 700 K. The dashed lines are 
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FIG. 4: (Color online) (a) Temperature dependence of the various energy modes of a suspended graphene sheet which is 
suspended along arm-chair direction, (b) Variation of molar heat capacity at constant volume for graphene subjected to 
compressive and tensile strains. 
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FIG. 5: (Color online) Temperature dependence of the bend- 
ing modulus of graphene (a) and variation of the nanotube 
curvature versus nanotube index, l.e m. 
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FIG. 6: (Color online) The absolute value of the square of the 
Fourier transform of atomic heights of C-atoms (l/iqp) versus 
the absolute value of the wave vectors of the graphene lattice. 



power law fits. Notice that a ^ —4 which clearly in- 
dicates that anharmonic effects are present in the used 
VFF potential. Moreover we see that a decreases with 
decreasing temperature (a =-3.0 for T=700K and a —- 
3.557 for T=200K) which hints that a more harmonic 
behavior is found at low temperature when using the 
VFF potential. The latter temperature dependence is 



in agreement with the REBO predictions [21|. Notice 
that the REBO is a bond-order interatomic potential. 
Note that the peaks in Fig. IH] are related to first Bragg- 
peak, 47r/3ao — 2. 94 A due to the discreteness of the 
graphene lattice. Notice that the modulation am- 
plitude in Fig. 6 for T=200 K is about O.sA and 
for T=700 K is about 0.7A which are temper- 
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ature and size dependent quantities, the larger 
the size the larger the amplitudes (here graphene 
has dimension 221 x 230 A-^). Here, we did not 
study the eflfect of size and refer the reader to 
Refs. \^, 10] where such a study can be found. 



VIII. CONCLUSIONS 

In this study, we showed that the recently proposed 
valence force field model (VFF) |14| for graphene en- 
ables to compare the contribution of the different energy 
terms when straining graphene. In a stretching experi- 
ment the main energy contributions are due to stretch- 
ing and bending terms while for compressive strains also 
other terms such as out-of-plane and tt — tt interaction 
terms play an important role. We found that using such 
a classical approach gives accurate values for the Young's 



modulus at room temperature which are found to be as 
accurate as those using ab-initio methods. The calcu- 
lated Young's modulus is close to the experimental re- 
sult. The total energy is quadratic in e for strains smaller 
than |2%|. The current VFF model predicts a temper- 
ature independence bending modulus. The temperature 
dependence of the total strain energy yields an accept- 
able value for the molar heat capacity of graphene which 
is too a large extend independent of the applied strain. 
The Gruneisen parameter is found to be positive and 
about 0.64. 
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